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Abstract 

We perform molecular dynamics and Monte Carlo simulations of two-dimensional melting with 
dipole-dipole interactions. Both static and dynamic behaviors are examined. In the isotropic liquid 
phase, the bond orientational correlation length £q and susceptibility X6 are measured, and the data 
are fitted to the theoretical ansatz. An algebraic decay is detected for both spatial and temporal 
bond orientational correlation functions in an intermediate temperature regime, and it provides 
an explicit evidence for the existence of the hexatic phase. From the finite-size scaling analysis 
of the global bond orientational order parameter, the disclination unbinding temperature Tj is 
estimated. In addition, from dynamic Monte Carlo simulations of the positional order parameter, 
we extract the critical exponents at the dislocation unbinding temperature T m . All the results are 
in agreement with those from experiments and support the KTHNY theory. 

PACS numbers: 64.70.Dv, 64.60.Ht 



1 



I. INTRODUCTION 



Two-dimensional melting has been intensively studied in the past years, but it is still 
not completely understood Melting in two dimensions is quite different from 

its counterpart in three dimensions, for a true long-range positional order doesn't exist 
in two-dimensional systems. The absence of a conventional long-range order at non-zero 
temperature was first pointed out by Mermin and Wagner J4]. Nevertheless, another long- 
range order, which is called the bond orientational order, can be observed in the solid phase 

0. 

There exist several possible theoretical descriptions of melting in two-dimensional sys- 



7J, |8| , predicts that 



terns. The KTHNY theory, developed by Halperin, Nelson and Young 
a third phase, the so-called hexatic phase, may exist between solid and liquid states in a 
portion of the phase diagram. The system first melts from the solid state to the hexatic state 
due to the unbinding of dislocation at a temperature T m , and then melts from the hexatic 
state to the liquid state at the disclination unbinding temperature Tj. Both transitions are 
Kosterlitz-Thouless phase transitions 3|. Naturally, the KTHNY theory only describes a 
possible scenario. It is also possible that anyone or both of the continuous transitions are of 
first order, and even that there is a direct first-order transition from the solid state to the 
isotropic liquid state. 

yen though quite some experiments supported the KTHNY theory 1(3, 
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171 ]. most early works of computer simulations on two-dimensional melting favored a 



first-order phase transition, and the hexatic phase was not observed. For example, for the 
systems with dipole-dipole interactions, Kalia and Vashishta |l8| observed a superheating 
and supercooling, as well as a latent heat in two-dimensional melting, and concluded that 
the phase transition is of first order. Later, Bedanov, Gadiyak and Lozovik [19( found that 
both the positional and bond orientational order vanished simultaneously at the melting 
point, and the hexatic phase didn't exist. Similar results have been found for other systems 
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23]. Even for the simp 



about the melting mechanism 
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est system, the hard disk model, there was no consensus 
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27]. 



In 1996, Bagchi et al. 28j carried out a finite-size scaling analysis of the bond orientational 
order parameter in a system interacting via a repulsive 1 / r 12 potential, and found that the 
results were in agreement with the KTHNY theory, even though no conclusive evidence 
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for the hexatic phase was observed. Later, extensive Monte Carlo simulations of the hard 



disk model were performed by Jaster 



29 



30J. Numerical behaviors of the susceptibility, 



spatial bond orientational correlation length and pressure, support the KTHNY theory. 
3ut the algebraic decay of the bond orientational correlation function was still not shown 



31|| . Recently, Monte Carlo simulations of a two-dimensional electron system with a 1/r 



interacting potential have been performed by He et al. 32|. An algebraic decay of the 
bond orientational correlation function is observed, and it explicitly reveals the existence 
of the hexatic phase. In principle, however, the finite-size effect and coexistence of liquid 
and solid phases may also lead to such an algebraic decay. One needs to carefully rule 
out these possibilities. On the other hand, in all these numerical simulations of the bond 
orientational order, the static behavior of the melting is mainly concerned, and the dynamics 
is not touched so much. 

Recently, more progress in computer simulations has been made in understanding two- 



33 



341 and exter- 



dimensional melting, for example, on the roles of the polydispersity 
nal fields 35], and on the structural change during the melting 36]. Especially, some 
experiments show much interest in a two-dimensional system with dipole-dipole interactions 



13, 



37] . The algebraic decay of the spatial and temporal correlation functions are observed 



and the dynamic behavior is found to be very relevant for two-dimensional melting. From the 
view of numerical simulations, the two-dimensional system with dipole-dipole interactions is 
not much understood. The purpose of this article is to present systematic computer simula- 
tions of two-dimensional melting with dipole-dipole interactions. Main results are obtained 
with molecular dynamics simulations, and Monte Carlo simulations are also performed in 
some cases and for confirmation. Both static and dynamic behavior will be examined, and an 
emphasis is given to the algebraic decay of both the spatial and temporal bond orientational 
correlation functions in the hexatic phase. 

The article is organized as follows. In Sec. II, the model and numerical methods will 
be described. In Sec. Ill, numerical results will be presented for both static and dynamic 
behavior. Finally it comes the conclusion. 
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II. MODEL AND METHOD 



A. The model 



In this article, we consider a two-dimensional dipolar system whose Hamiltonian can be 
written as 

i 2 i=1 \rij\ 3 \ r ij\ 5 

where pi, fii and frti are the momentum, mass and magnetization of the zth dipole respec- 



tive 



13, 



y, and N is the total number of particles. In order to mimic the experiments in Refs. 



371 ] and to simplify the problem, we assume the dipoles are aligned perpendicular to the 



surface. Thus, Eq. ([I]) can be reduced to 



N 2 i N N 

^=e^ee9S- (2) 

i * i=i i ' ij i 



For convenience in numerical simulations, we rewrite Eq. (T5]) as 



N v 2 1 N N a 

^ = E| + iEE^) 3 , (3) 

where we have assumed the mass and magnitude of the magnetization of the dipoles are 
identical. For simplification, the reduced units are adopted, in which the parameters e 
and a, Boltzmann constant and mass [i of the dipoles are set to 1. The thermodynamic 
observables are determined only by a dimensionless constant T = ecr 3 (7m) 3 / 2 jkT [38(] , where 
n = N/V is the 2D volume fraction of the dipoles. 

The reasons we choose this model are: (i) this model lacks extensive numerical study, and 
the existing results do not favor the KTHNY theory; (ii) there are unambiguous experimental 
results of such system [17j, |37|, to which we may compare our results. 



B. Numerical methods 

In our simulations, particles are put in a rectangular box with a size ratio 2 : the 
density of the particles is fixed to be l/(2y / 3), and the number of the particles is taken to 
be from 1024 to 32768. The linear size L of the system is related to the total number iV 
of particles by L = Periodic boundary conditions are used in simulations, and the 

dipole-dipole potential is truncated at 10. In two dimensions, such a truncation is reasonable. 
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Actually, the correction of the potential to the truncation is u c = J 27r J^°° g(r)e(^) 3 rdrd6; 
assuming g(r) = 1, it leads to u c = 2irea 3 /r c , which decays to zero with r c . In fact, the main 
parts of the simulations are carried out in the hexatic and liquid phases where g(r) quickly 
stabilizes at a constant which is smaller than 1 (see Fig. H] (b)). 

In order to confirm the truncating procedure, we have performed the simulations at dif- 
ferent truncating distances, and find that the difference is neg ligibly small. In addition, 
extra simulations using the Ewald Summation technique 39_|, |40| , which is known for trans- 
ferring long-range interactions to short-range ones, are also performed to further justify our 
truncation. Within statistical errors, the results for the global bond orientational order pa- 
rameter \l/ 6 and susceptibility Xe obtained with different truncating distances and the Ewald 
Summation are in good agreement with each other. Relevant data with L = 64 are compiled 
in tabled for comparison. Additional simulations with L = 128, and measurements of the 
pair distribution function g(r) also confirm the reliability of the truncation. 

In this paper, most simulations are performed with molecular dynamics. All results 
are obtained at a constant temperature with the NVT ensemble based on the Nose-Hoover 
Chain thermostat [44, 42]. The equation of the motion is solved via the five-point Nordsieck- 



Gear predictor-corrector method. The time step At in all the simulations is set to 0.01. A 
shift of the conserved total energy is within 0.0001%. 

The initial configurations in our simulations consist of particles uniformly distributed 
over the system box on a triangular lattice. Before collecting data for the measurements of 
physical observables, the system is carefully equilibrated, especially in the critical regime. 
We monitor the global bond orientational order parameter, and begin our measurement 
after this order parameter reaches a steady value. For the larger system (N = 16384), for 
example, it takes 5 x 10 5 time steps to thermalize the system. Only the configurations in 
equilibrium are used for the measurements of observables, extending over 9 x 10 6 time steps. 
In order to obtain independent configurations, the autocorrelation function of the global 
bond orientational order parameter is measured, and the correlation time is estimated to be 
r m 2400 time steps in the critical regime. Then the measurement is performed every 2500 
time steps. 

In order to confirm our molecular dynamics simulations, standard Monte Carlo simu- 
lations are additionally performed. For example, The bond orientational correlation func- 
tions from both molecular dynamics simulations and Monte Carlo simulations are shown in 
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Fig. [T] (a). Both methods provide consistent curves, and it shows that our molecular dy- 
namics simulations indeed generate proper ensemble distributions. Furthermore, dynamic 
Monte Carlo simulations are carried out to extract the critical exponents for the positional 
order parameter at the dislocation unbinding temperature T m . 



C. Observables 



The bond orientational symmetry of a solid can be described by the six-fold global bond 
orientational order parameter \l/ 6 defined as 

1 N 

*« = <It;£iM>, (4) 
iv k=l 

where N is the total number of the particles, (• • ■) denotes the ensemble average or the time 
average in molecular dynamics simulations and Monte Carlo simulations, and the ipg^ is the 
local bond orientational order parameter defined as 

= ^rE ex P( i6 M- ( 5 ) 

k j 

Here the sum j is over the neighbors of the particle k, and 8kj is the angle between r^} (the 
relative position vector of the particle k and j) and an arbitrarily fixed reference axis. Neigh- 



bors are obtained with the Voronoi polygon 43] . The susceptibility of the bond orientational 
order is defined as 

X* = N(*l). (6) 

The hexatic phase is characterized by an algebraic decay of the bond orientational cor- 
relation function defined as 

9e(ri -r$) = (^(^1)^(^2))- (7) 

In order to obtain an accurate value of the bond correlation length, we smooth the bond 
orientational correlation function following Ref. 30]. We divide the volume of the system 
into stripes with a width of Ax and measure the bond orientational correlation between 
different stripes. 

1 rL px+Ax/2 _^ 1 pL /•Ax/2 __x 

%{x) = ((— - dy' d*V 6 , fe (?)r x (— - / dW d*V 6 , fc (7))>, (8) 

J\{X) JO Jx-Ax/2 N{0) JO J-Ax/2 
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where 

i-L rx+Ax/2 _j> 

N(x) = dy' dx'p(r), (9) 

JO Jx-Ax/2 
N 

p ^ = J2S{-t-rt), (10) 

i=i 

and L is the linear size of the system in the y direction. The temporal bond orientational 
correlation function characterizes the time correlation of the bond orientational order pa- 
rameter, and is defined as 

9e(t) = &ik(t°)Mto + t)), (11) 

where ip6,k{to) and i^>e,k{to + 1) are the local bond orientational order parameters measure at 
the time to and to + t respectively, and the average is over to i n equilibrium. In the hexatic 
phase, ge(t) also decays by a power law [371 ]. 

The positional symmetry of solid can be described by a positional order parameter defined 

as 

1 N -> 

% 0S = {^J2exp(ia-r])), (12) 
j=i 

where G is a reciprocal-lattice vector which gives the first Bragg peak. In practice, we 
average the order parameter over the six reciprocal vectors which correspond to the six 
vectors connecting the six neighbors from the lattice site j. The positional correlation 
function is defined as 

9G{\t ~ P\) = (exp0 ■ ( r> - ?))). (13) 

In the hexatic phase, the positional correlation function decays exponentially. Finally, the 
pair distribution function is defined as 

5( r ) = ^E^-4 ( 14 ) 



N 2 . 



where V is the volume of the system. 



III. COMPUTER SIMULATION 



In this article, we perform extensive simulations of two-dimensional melting in the NVT 
ensemble with system sizes up to 32786 atoms, and find a strong evidence for the existence 
of the hexatic phase in the dipole-dipole interacting system. We first measure the spatial 
bond orientational correlation function and susceptibility in the isotropic liquid phase and 
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compare the results with the predictions of the KTHNY theory. This gives us estimates 
of the isotropic- anisotropic transition temperature Tj. With this critical temperature in 
hand, we further scan the parameter space, and observe an algebraic decay of the spatial 
bond orientational correlation. We also measure the temporal bond orientational correlation 
function, and its behavior is in good agreement with the KTHNY theory. In order to rule 
out a possible coexistence phase and the finite-size effect, we perform a homogeneous test 
and finite-size scaling analysis of the bond orientational order parameter. The result is 
compatible with previous measurements. At last, with Monte Carlo methods we simulate 
the short-time dynamics of the positional order and estimate the exponent f] m , and the value 
is also in agreement with the theoretical prediction. All our results are compatible with the 
experiments and KTHNY theory, and the hexatic phase is explicitly observed. 



A. Bond orientational order 



The bond orientational order parameter ^/q offers a direct description of the bond orienta- 
tional order 23j . Assuming Tj is the transition temperature of the bond orientational order 



and T m is the transition temperature of the positional order, the bond orientational order 
parameter should vanish for T > Tj. and take a finite value less than 1 for T < T m . The 
behavior of at the temperatures between Tj and T m depends on the underlying melting 
scenario. If the transition at T is of first order, increases linearly from Tj to T m . If 
the melting scenario is of KTHNY, i.e., the transition at Tj is a Kosterlitz-Thouless phase 
transtion, \l/ 6 then vanishes throughout the hexatic phase for there doesn't exist a true long- 
range bond orientational order. However, the finite-size effect in the simulations blurs this 
distinction and prevents us drawing a clear conclusion. Nevertheless, the measurement of 
the bond orientational order parameter \l/ 6 does give us an estimated value of Tj « 0.01250. 
In the Fig. [1] (b), ^§ versus T is shown. 

To further understand the phase transition at Tj, we measure the correlation length and 
susceptibility of the bond orientational order parameter in the isotropic liquid phase for 
different temperature T. For the measurements are carried out in the isotropic liquid phase, 
the spatial bond orientational correlation function is independent of the spatial directions. 
We extract the correlation length £ from the exponential decay of the bond orientational 



S 



correlation function smoothed with the technique described in Eq. (jSJ), 

g 6 (x) ~ exp(-z/0- (15) 

Subsequently, we compare our results with the predictions of the KTHNY theory, i.e., an 
exponential singularity of the correlation length and susceptibility, 

e 6 (r) = a € exp(6 ? r- 1 / 2 ), (16) 

Xe(r) = a x exp(6 x r" 1/2 ), (17) 

as r = T — Tj — > + . In Fig. [21 the numerical data are fitted to the above exponential 
forms. The best fit of the correlation length and susceptibility gives Tj = 0.01237(16) and 
Tj = 0.01243(4) respectively. These two values are in agreement with each other within 
statistical errors, and are also consistent with the previous estimated value of Tj from the 
global bond orientational order parameter. Therefore, our results support the KTHNY 
prediction, even though the statistical errors of the correlation length and susceptibility are 
relatively large. 

B. The hexatic phase 

According to the KTHNY theory, the hexatic phase is characterized by an algebraic decay 
of the bond orientational correlation function and an exponential decay of the positional 
correlation function. Therefore, we scan the parameter space between Tj and T m , and 
measure the bond orientational correlation function and positional correlation function. The 
bond orientational correlation function is shown in Fig. [3] (a). A clear evidence for the 
existence of the hexatic phase is observed. 

i) In the solid phase (T = 0.0119), the correlation function rapidly stabilizes to a constant, 
and it indicates that there is a true long-range order of the bond orientational symmetry. 

ii) In the hexatic phase (T = 0.01252 and 0.01253), the correlation function shows an 
algebraic decay, 

3 6 (x)~x-" 6 , (18) 

and it indicates that there is a quasi-long-range order of the bond orientational symmetry. 

iii) In the liquid phase (T = 0.0127), the correlation function decays exponentially, and 
it indicates an isotropic state. 
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After smoothing the correlation function at T = 0.01253 , we obtain an exponent rje = 
0.252(6) from the slope of the curve, and it is close to the value rje = 0.25 at the transition 
temperature Tj predicted by the KTHNY theory. If we assume T = 0.01253 is just the 
transition temperature Tj, it is quantitatively in agreement with the previous measurements 
in the preceding subsection. Nevertheless, to obtain a more accurate value of Tj , we need 
to consider more carefully the finite-size effect. In the next subsection, we will locate the 
transition temperature from the finite-size scaling. 

In Fig. HI the positional correlation function and pair distribution function are shown 
respectively. One may observe two different behaviors. 

i) In the solid phase (T = 0.0119), the positional correlation function shows an algebraic 
decay, indicating a quasi-long-range positional order in a two-dimensional solid, while the 
oscillation in the pair distribution function persists over the entire range. 

ii) In the hexatic phase (Tj = 0.01253), the positional correlation function decays quickly 
to zero, indicating that there exists no positional order in the hexatic phase, and the os- 
cillation in the pair distribution function dies off rapidly. The behaviors of the positional 
correlation function and pair distribution function in the liquid phase are qualitatively the 
same as in the hexatic phase. 



The recent experiment reported in Ref. 37] shows that the dynamic behavior is also 



very important in understanding the melting mechanism in two dimensions. According to 
the KTHNY theory, in the solid phase the temporal bond orientation correlation function 
will rapidly stabilize to a constant, in the hexatic phase it shows an algebraic decay with an 
exponent equal to fje/2 , 

gs{t)~r^\ (19) 



and in the liquid phase the temporal correlation function decays exponentially [441 ] . Such a 
behavior is indeed observed in experiments, and it provides a strong evidence for the exis- 
tence of the hexatic phase. To our knowledge, such measurements have not been performed 
in numerical simulations. 

In order to deepen our understanding of two-dimensional melting and further confirm 
our observation in numerical simulations of static properties, the temporal bond orienta- 
tion correlation function is measured in our molecular dynamics simulations. The result is 
shown in Fig. [3](b). Obviously, as the temperature changes from T = 0.0119 to 0.01253, 
then to 0.0131, the temporal bond orientation correlation function follows the prediction of 
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the KTHNY theory, and are well consistent with the experimental observation 37|]. The 
exponent measured from the slope of the curve at T = 0.01253 is 0.0843, somewhat smaller 
than the theoretical prediction 0.125 at Tj. This probably suggests that the anisotropic- 
isotropic transition temperature 7] should be still slightly above the value 0.01253, and our 
measurements of the spatial and temporal bond orientational correlation functions may still 
carry certain finite-size effects. 



C. Finite-size scaling analysis 

Our measurements of the spatial and temporal bond correlation functions provide us an 
explicit evidence for the existence of the hexatic phase in two-dimensional melting. However, 
it is difficult to extract an accurate transition temperature Tj from the correlations functions. 
One may directly measure the correlation length in the isotropic liquid phase and then fit 
the data to the ansatz in Eq. ( TT6l) and obtain the transition temperature Tj. But this 
approach also suffers from the difficulty that one can not obtain the correlation length 
at the temperatures very close to Tj 45, [46] . Meanwhile, due to the finite-size effect, 



distinguishing between an algebraic and an exponential decay might be problematic if the 
correlation length is finite but much larger than the system size. Therefore, to extract a more 
accurate disclination unbinding temperature Tj and to confirm the previous observation of 
the hexatic phase, we perform a finite-size scaling analysis of the bond orientational order 
parameter. 

From the finite-size scaling form, the second moment of the bond orientational order 
parameter can be written as 

(*g>~L-*/(L/e 8 ), (20) 

where L = 2\^N is the linear size of the system and is the bond correlation length. Since 
^6 is divergent in the hexatic phase, (^g) thus shows a power-law behavior with respect to 
L in the hexatic phase. In the liquid phase, the power-law behavior will be modified by the 
scaling function f(L/^ e ). 

We measure the second moment of global bond orientational order parameter with system 
size L = 64, 128, 256 at T = 0.01257 and perform finite-size scaling analysis mentioned 
above. The open circles shown in Fig. [5] are the results. In order to save computation time, 



we use the subsystem method introduced by the authors of Refs. [23 



Here, a brief 
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comment about the above non-standard finite-size scaling analysis is needed. In principle, 
the subsystem procedure still carries a second-order finite-size effect induced by the finite 
bulk system size L. But this second-order finite-size effect is negligibly small in practical 
simulations [23j, and the procedure has been proved to be reliable and may reduce computer 



times [28J. We also carried out the finite-size scaling analysis using subsystem method at 
T = 0.01257 to further justify this procedure, the result is shown in Fig. [5j It is easy to 
observe that within statistical errors, the data with periodic boundary conditions and from 
subsystems are well consistent. With the subsystem method, we measure (\I/q) at different 
temperatures with a bulk linear size L = 256 or 512, and a total number of particles ranging 
from N = 16384 to 32768. Then the system is divided into small subsystems with a linear 



size Lg 23j and the bond orientational order parameter of each subsystem is measured. The 
result is shown in Fig. [5j 

To locate Tj, we assume T]q = 1/4. In other words, we search for a temperature which 
yields t/q = 1/4, and then assign this temperature to be Tj. The requirement of tjq = 1/4 



yields the upper limit of Tj 44]. Combining the results obtained in the preceding subsections, 
we conclude 0.01253 < T < 0.01257. To compare our results with those in literatures, we 
convert T to the dimensionless parameter Tj, and obtain 68.707 < Tj < 68.927. It improves 
the values Tj = 62 ± 3 with a small system_ N = 256 [18j and T = 67.750 with a relatively 
larger system N = 961 j47j. In Refs. 



18 



471 ] , the phase transition is supposed to be of first 
order, and the values of Tj are obtained from the hysteresis in the temperature dependence 
of energy, the existence of latent heat and the thermodynamic nucleation of the solid from 
the supercooled liquid. Our estimate of Tj is based on the KTHNY theory, and is much less 
affected by the finite-size effect. 

D. Ruling out the coexistence phase 

In principle, the NVT molecular dynamics simulation can not obviate the coexistence 
phase, and the superposition of the solid and liquid phases may also produce the hexatic-like 
behavior. In order to exclude this possibility, we apply the homogeneous test. We divide 



the system into small subsystems and compute the susceptibility x& f° r & U subsystems 21]. 
If the system exhibits an inhomogeneous two-phase coexistence, the probability distribution 
of xg & t a sufficiently small length scale could be modeled by a curve with two peaks, which 
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reflects a combination of solid and fluid distributions. On the other hand, if the system is 
homogeneous, varying the size of the subsystems should not lead to any qualitative change 
in the probability distribution of Xsi i- e -> the curve should always remain with a single peak. 

We have measured the possibility distribution of Xe m the hexatic phase at T = 
0.01257, 0.01253, 0.01252 , and in order to compare the result with that in the homogeneous 
phase, we also perform a simulation at an extra temperature T = 0.0100 corresponding to 
the cool solid phase. No qualitatively change is found at these temperatures. This test rules 
out the existence of a coexistence phase, and confirms the observation of the hexatic phase 
in the previous subsections. The result at T = 0.01252 is shown in Fig. El It is clearly seen 
that varying the system size doesn't change the shape of distributions, but only shifts the 
peak of the probability distribution. 



E. Dynamic Monte Carlo simulations 



In the last decade, it has been discovered that already in a macroscopic short-time regime 



emerges the universal scaling behavior 



48 



49, 



50 



51 



52j. Measurements now are carried out 



at the early stage of the time evolution, therefore one does not suffer from critical showing 
down. The dynamic scaling form of the second moment of the positional order parameter 
below the dislocation unbinding transition temperature T m is 



*Ut,L)=b-^ 2 (b-%b- 1 L) 



(21) 



where t is the evolution time, z is the dynamic critical exponent, and b is an arbitrary 
rescaling factor. For a sufficient large L, this dynamic scaling form is reduced to 



(22) 



From a finite-size scaling analysis of the time-dependent Binder cumulant 48] 



Upos (^) 



\]> 4 



( ^ pos ) 



one obtains 



U pos {t)~t d ' z /L a 



(23) 



(24) 



where d is the dimension of the system. The dynamic critical exponent z can be estimated 
from Eq. fl24|) . and with z in hand, the static exponent rj m can be obtained from Eq. fl22|) . 
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Now, we turn to locate the transition temperature T m . As the temperature increases, 
T m is characterized by the dislocation unbinding which breaks the quasi-long-range posi- 
tional symmetry. Therefore, one may measure the correlation function of the positional 
order parameter to estimate T m , for the positional correlation become short-range at T m . 
Nevertheless, this method suffers from the difficulty that one needs to do simulations in 
the critical region. Even for the hard disk model, in which the thermodynamic quantities 
are only determined by the density p of disks, it is still not easy to locate p m accurately. 
p m « 0.933 is reported in Ref. 50], while p m = 0.910(2) is given in Ref. [521 ] . 

Alternatively, a dynamic technique for locating T m is applied in the experiments reported 



371 ] and adopt the dynamic 



in Ref. [37j. In our computer simulations, we follow Ref. 
criterion for T m , since the method is relatively simple, and may provide direct comparison 
with experiments. We first introduce the 2D Lindemann parameter 
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53] 



l?(~r* + Oq) — Ut^)] 2 ) x Tin, 



(25) 



where aTj is the lattice spacing vector, ~r* is the positional vector, ~Tt is the displacement field 
and n is the 2D volume fraction of particles. Initially, the system is set on a perfect trian- 
gular lattice. In numerical simulations, we gradually warm up the system with the velocity 
rescaling procedure. At each T, the system is equilibrated to the thermal equilibrium. Then 
we measure the Lindemann parameter at different temperatures. In general, a sharp growth 
of 7 m indicates vanishing of the positional symmetry. Such a Lindemann criterion in 3-D 
systems is a well-established and justified procedure for locating the melting temperature 
T m , although it was unclear in two dimensions [4j. In 1985, Bedanov and Gadiyak improved 
the definition of the Lindemann parameter to the form in Eq. (|25l) and demonstrated in the 
simulations of electron and dipole systems that when 7 m rises up to a critical value 0.12, 
the melting takes place js^j]. At the melting point T m , which is more clearly identified by 
the sudden drop of the positional correlation length, a sharp growth of 7 m is induced by 
the leap of disclination number and self-diffusion constant. Therefore, this local quantity is 
relevant to the melting. Recently, Zahn et al. applied this criterion to their experiments 
0, 0], and the results are in agreement with those from numerical simulations 53]. The 
Lindemann parameter may provide at least a first estimate of the melting temperature T m . 
Due to its efficiency and simplicity, the Lindemann criterion has been applied to different 



systems for locating the melting point 



54 



55 



56] 
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Here we locate T m with the Lindemann criterion. Four runs are performed in order 
to estimate T m . One of them is shown in Fig. We estimate T m = 0.0120(2). The 
hexatic phase of the dipolar system lies in a range of the phase diagram, between 0.01253 < 
Tj < 0.01257 and T m = 0.0120(2). This is comparable with that of the hard disk model, 
Pi = 0.899(1) 0, and p m « 0.933 in Ref. Q while p m = 0.910(2) in Ref. 52]. T t and 
T m may overlap for small systems, and this is one reason why the hexatic phase was not 
observed in some previous studies. 

Now we perform dynamic Monte Carlo simulations at the transition temperature T m . 
The reason we perform Monte Carlo simulations is that the dynamic scaling forms in Eqs. 
(T2"2"j) and (1241) may not hold in the dynamic process of Nose-Hoover chain molecular dy- 
namics simulations. It seems that the Nose-Hoover Chain method is originally devised for 
equilibrium simulations and contains techniques violating the dynamic scaling behavior. 
In comparison to this, the dynamic scaling behavior in Monte Carlo simulations has been 
extensively justified. 

In Monte Carlo simulations, the system initially at an ordered state is released to the 
dynamic evolution with the Metropolis algorithm, and then the time-dependent ^ pos and 
U pos are measured. By fitting ^/ pos (t) and U pos (t) to Eqs. fl22l and fl24"|) . both the dynamic 
exponent z and static exponent r\ m can be determined. For comparison, we also perform the 
same simulations at another temperature T = 0.0115. The results are shown in Fig. [8] (a) 
and (b). 

From U poa (t) in Fig. M (b), we estimate z = 1.910(70), and from ^l os (t) in Fig. M (a), we 
measure r\ m jz = 0.143(5). Combining these results, we deduce r) m = 0.273(20). This value 
is also in agreement with the prediction (1/4 < r] m < 1/3 ) based on the KTHNY theory 

A- 



IV. CONCLUSION 



We present molecular dynamics and Monte Carlo simulations of two-dimensional melting 
with dipole-dipole interactions. An algebraic decay is observed for both the spatial and 
temporal bond orientational correlation functions in an intermediate temperature regime, 
and this serves as an explicit evidence for the existence of the hexatic phase. 

To obtain a relatively accurate disclination unbinding temperature Tj, we perform a finite- 
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size scaling analysis for the bond orientational order parameter. The result 0.01253 < T; < 
0.01257 improves the value from a direct fit of the correlation length to the exponential 
ansatz. In addition, by analyzing the probability distribution of the bond orientational 
susceptibility x 6 , a possible coexistence phase is ruled out. 

At last, from dynamic behavior of the Lindemann parameter, the dislocation unbinding 
transition temperature is estimated to be T m = 0.0120(2). We also perform dynamic Monte 
Carlo simulations of the positional order parameter and the time-dependent cumulant. From 
the power-law behavior of these quantities, we determine the exponents r\ m = 0.273(20) and 
z = 1.910(70) . 

In summary, a clear evidence for the existence of the hexatic phase is observed for two- 
dimensional melting with dipole-dipole interactions, and all the static and dynamic behaviors 
of the system are compatible with recent experiments and the KTHNY theory. 
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# 6 (T = 0.0150) 


^ 6 (T = 0.0125) 


X6 (T = 0.0150) 


Xe (T = 0.0125) 


r t = 10 


0.0842(25) 


0.684(4) 


9.14(49) 


479(6) 


r t = 20 


0.0849(35) 


0.680(4) 


9.34(38) 


475(5) 


Ewald Summation 


0.0859(37) 


0.680(1) 


9.50(84) 


474(2) 



TABLE I: The global bond orientational order parameter \I>6 and susceptibility xe measured by 
truncating the potential at r t = 10, r t = 20 and with the Ewald Summation to deal with the 
potential. The linear size is L = 64, and the temperature is T = 0.0150 and T = 0.0125. 
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FIG. 1: (a) ge(x) obtained with molecular dynamics (MD) and Monte Carlo (MC) simulations at 
T = 0.0150 plotted vs. £ on a semi-log scale. The smoothed curve is shifted upward for clarity. 
The smoothing technique is described in Sec. II. (b) ^6 plotted vs. T on a linear plot. The bond 
orientational order parameter increases abruptly around T = 0.0125. The line fitted to the circles 
is a guide to the eyes. 




FIG. 2: Bond orientational correlation length (full symbols) and susceptibility (open symbols) as 
a function of temperature. The curves show the best fits of Eqs. (|16p and (|17p according to the 
KTHNY theory. The fitted transition temperatures are T { = 0.01237(16) and 0.01243(4) for the 
correlation length and susceptibility respectively 
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FIG. 3: (a) The spatial bond orientational correlation function g® (x) plotted vs. x on a double 
decimal log scale. The temperature T = 0.0119 is just before melting, T = 0.0127 is typically in 
the liquid phase, and T = 0.01252 and 0.01253 are in the hexatic phase. The straight line with 
a slope of —1/4 is a guide to the eyes, (b) The temporal bond orientational correlation function 
<76 (t) plotted vs. t on a double-log scale. The temperature T = 0.0119 is just before melting, 
T = 0.0131 is typically in the liquid phase, and T = 0.01253 is in the hexatic phase. ge(t) at 
another T = 0.01257, which is slight above the estimated Tj, is also shown. Lines fitted to the data 
are to guide the eyes. 
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(a) (b) 

FIG. 4: (a) The positional correlation function gc{x) at T = 0.01253 in the hexatic phase (lower 
curve) and 0.0119 in the solid phase (upper curve) plotted vs. x on a linear scale, (b) The pair 
distribution function g(x) plotted vs. x on a linear scale. The upper two curve are shifted upward 
for clarity. The curves at T = 0.0119, T = 0.01253 and T = 0.0127 show features in the solid, 
hexatic and liquid phases. 




FIG. 5: The finite size scaling analysis of ^q. L = 256 or 512 is the bulk linear size and L s is the 
size of the subsystem. The dotted line with a slope of —1/4 is a guide to the eyes. Open circles are 
the results from independent simulations with periodic boundary conditions at L = 64, 128, 256. 
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FIG. 6: The probability distribution of %6 m the hexatic phase at T = 0.01252. The symbols in 
the figure indicate the mean numbers of particles in different subsystem. 
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FIG. 7: The Lindemann parameter 7 m vs. T. The critical temperature T m = 0.0120 is visualized 
by the vertical dotted line. 
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FIG. 8: (a) ^fp 0S (t) plotted vs. t on a double-log scale. The lower curve is shifted downward for 
clarity, (b) U pos {t) plotted vs. t on a double-log scale. The upper curve is shifted upward for 
clarity. 
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